# Install the libraries (if using binder)
#!pip install numpy
#!pip install pandas
#!pip install matplotlib
#!pip install seaborn
#!pip install pylab
#!pip install plotly
#!pip install geopy

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pylab as plot
import plotly.graph_objects as go
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="my_user_agent")

Reading the csv file and get it into a dataframe format¶

df = pd.read_csv("pollution_us_2000_2016.csv")

print(f'The dataframe has {len(df)} rows and {df.shape[1]} columns')

The dataframe has 1746661 rows and 29 columns

Based on the above numbers, one concludes that the dataframe corresponds to very large dataset (Big Data). Therefore, it is not advisable to view the entire dataframe as that will be a super memory expensive task. One needs to perform the analytics in a smarter way so as to the get the desired results while making sure not to put too much load on the memory.

all_states = df['State'].unique()
num_states = len(all_states)

print(f'There are {num_states} states in the dataset.')

There are 47 states in the dataset.

print(f'The columns in the dataframe are given by \n{df.columns}')

The columns in the dataframe are given by 
Index(['Unnamed: 0', 'State Code', 'County Code', 'Site Num', 'Address',
       'State', 'County', 'City', 'Date Local', 'NO2 Units', 'NO2 Mean',
       'NO2 1st Max Value', 'NO2 1st Max Hour', 'NO2 AQI', 'O3 Units',
       'O3 Mean', 'O3 1st Max Value', 'O3 1st Max Hour', 'O3 AQI', 'SO2 Units',
       'SO2 Mean', 'SO2 1st Max Value', 'SO2 1st Max Hour', 'SO2 AQI',
       'CO Units', 'CO Mean', 'CO 1st Max Value', 'CO 1st Max Hour', 'CO AQI'],
      dtype='object')

One would now like to get rid of the columns that are not going to be useful for our analytics.

# columns to be dropped
drop_cols = ['Unnamed: 0', 'State Code', 'County Code', 'Site Num', 'Address', 'NO2 1st Max Value', 'NO2 1st Max Hour', 'NO2 AQI', 'O3 1st Max Value', 'O3 1st Max Hour', 'O3 AQI', 'SO2 1st Max Value', 'SO2 1st Max Hour', 'SO2 AQI', 'CO 1st Max Value', 'CO 1st Max Hour', 'CO AQI']

df.drop(drop_cols, axis=1, inplace=True)

print(f'The dataframe now contains {len(df)} rows and {df.shape[1]} columns')

The dataframe now contains 1746661 rows and 12 columns

To get some idea about the entries in the dataframe, instad of looking at the entire dataframe, one peeks at only few rows.

df.head()

Conversion of Units to achieve uniformity¶

By peeking at the above tiny dataframe, one sees that the units are not uniform for all the gases. Before converting all the units, one must make sure if all the entries in particular units column are the same, i.e., for a gas x, one would like to make sure if it is always measured in the same units.}

def same_units(df, col):
    """ Checks if all the entries of a column in the dataframe are same. In the context of the project,
    it checks if all unit columns have same units.
    
    Args:
        df (dataframe): the dataframe whose column is to be inspected.
        col (str): the column name whose entries are to to be inspected.
    
    Returns:
        True if all the entries are same, else False.
    """
    if len(df[col].unique()) == 1:
        return True
    return False

def get_units(df, col):
    """ 
    Returns all the units present in the unit column of the dataframe.
    """
    return df[col].unique()

same_units(df, 'NO2 Units')

True

get_units(df, 'NO2 Units')

array(['Parts per billion'], dtype=object)

same_units(df, 'SO2 Units')

True

get_units(df, 'SO2 Units')

array(['Parts per billion'], dtype=object)

same_units(df, 'O3 Units')

True

get_units(df, 'O3 Units')

array(['Parts per million'], dtype=object)

same_units(df, 'CO Units')

True

get_units(df, 'CO Units')

array(['Parts per million'], dtype=object)

Note that, $NO_2$ and $SO_2$ are always measured in parts per billion (ppb), while $O_3$ and $CO$ are always measured in parts per million (ppm).

A standard unit in the field of science is parts per million (ppm), hence one would like to convert $NO_2$ and $SO_2$ into parts per million (ppm).

df['NO2 Mean'] = df['NO2 Mean']/1000
df['SO2 Mean'] = df['SO2 Mean']/1000

The columns describing the units are now irrelevant, hence, one would like to discard all of these columns to get even more finer dataframe.

unit_cols = ['NO2 Units', 'SO2 Units', 'O3 Units', 'CO Units']
df.drop(unit_cols, axis=1, inplace=True)

Before, moving ahead, a good practice is to check if there are any NaN, i.e., missing values in the dataframe.

df.isna().any()

State         False
County        False
City          False
Date Local    False
NO2 Mean      False
O3 Mean       False
SO2 Mean      False
CO Mean       False
dtype: bool

Therefore, none of the above columns contain any missing value. Great!!

Rename the 'Date Local' column and make it the index of the dataframe.¶

To see the trends with time, it is a good practice to have the date column as the index of the dataframe.

df.rename(columns={'Date Local':'Date'}, inplace=True)

# first check the type of the Date column
df['Date'].dtypes

dtype('O')

One would like to convert this into a 'DateTime Object' since that makes it easier to parse dates and do analysis on it.

df['Date'] = pd.to_datetime(df['Date'])

# make 'Date' as the index of the dataframe
df.set_index('Date', inplace=True)

# let us take a peek at the dataframe
df.tail()

Insert Total Mean Column¶

In order to quantify $NO_2$, $O_3$, $SO_2$ and $CO$ as one single entity for each state, one would like to get a Total Mean column, that allows us to do exactly that. The entries of this columns will act as a a measure for each state and will help us see how each of these states have performed through time.

df['Total Mean'] = df['NO2 Mean'] + df['SO2 Mean'] + df['CO Mean'] + df['O3 Mean']
df.head()

Get individual Dataframe for each of the Gases¶

The goal now is to create individual dataframes for each of the gases, with columns representing each of the states. One would also like to discard the states that have missing value for any of the years.

Note that, we the original dataset has 47 states (as computed before). It might be the case, that not all of them have values recorded for all the yeras from 2000 to 2016. We would like to get rid of them.

def get_df(gas_name):
    """
    Computes dataframe for a particular gas with rows as years from 2000 to 2016 and columns as states.
    
    Args:
        gas_name (str): can be any of the NO2 Mean, SO2 Mean, O3 Mean, CO Mean or Total Mean.
    
    Returns:
        result (dataframe): rows are years, columns are states and an entry corresponds to average of all the 
                            emissions of the gas for the corresponding year and state. 
    """
    
    global df, all_states
    
    df_year_list = [df.loc[str(i)] for i in range(2000, 2017)]
    
    # populates dictionary with states as the keys and list of values
    dict_of_lists = {state:[] for state in all_states}
    
    for df_year in df_year_list:
        for state in all_states:
            # extract dataframe for that particular state
            state_vals = df_year[df_year['State'] == state]
            # get the average of all the values for the particular state and df_year
            dict_of_lists[state].append(state_vals[gas_name].mean())
    
    # convert the dictionary to a dataframe
    result = pd.DataFrame(dict_of_lists)
    
    # change the indices so that they reflect years
    index_list = [item for item in range(2000, 2017)]
    result.index = index_list
    
    # get only the columns for which one has non-NaN values.
    result.dropna(axis = 1, inplace=True)
    return result

# dataframe for no2
no2_df = get_df('NO2 Mean')

# dataframe for so2
so2_df = get_df('SO2 Mean')

# dataframe for co
co_df = get_df('CO Mean')

# dataframe for o3
o3_df = get_df('O3 Mean')

# dataframe for 'Total Mean'
total_df = get_df('Total Mean')

Let us check one of the dataframes to see how it looks:

no2_df

print(f'There are {len(no2_df.columns)} states in the datframe corresponding to the nitrogen dioxide emissions.')

There are 14 states in the datframe corresponding to the nitrogen dioxide emissions.

Naturally, one would like to know if all the dataframes have same states.

print(total_df.columns == no2_df.columns)
print(no2_df.columns == o3_df.columns)
print(o3_df.columns == so2_df.columns)
print(so2_df.columns == co_df.columns)

[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True]
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True]
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True]
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True]

Therefore, only these 14 states have their data recorded from 2000 to 2016. Rest other states have missing data, hence are discareded in our analysis.

Line Plots¶

Let us plot individual line plots for each of the gases to infer a visualization trend of how each of the 14 states are performed from 2000 to 2016 for each of the gases.

params = {'legend.fontsize': 23,
          'legend.handlelength': 2,
         'font.size': 27}
plot.rcParams.update(params)

# Setting the style of the Seaborn line Plots
sns.set_style("whitegrid")
plot_colors = ['#ff0000', '#940606', '#5c1414', '#ff8c00', '#965709', '#b59509', '#91790c', '#84910c', '#c1db00', '#37db00', '#216e08', '#00ffd9', '#033b6b', '#0089ff']
plot_line_styles = ['-', '--','-.', '-.', '--', '-', '--','-.', '-.', '--', '--','-.', '-.', '-']

# SO2
so2_df.plot(figsize=(40,25), color = plot_colors, linewidth=3.5, style = plot_line_styles)
plt.title("Sulphur Dioxide Trends from 2000 to 2016")
plt.xlabel("Year")
plt.ylabel("SO2 (ppm)")

Text(0, 0.5, 'SO2 (ppm)')

# NO2
no2_df.plot(figsize=(40,25), color = plot_colors, linewidth=3.5, style = plot_line_styles)
plt.title("Nitrogen Dioxide Trends from 2000 to 2016")
plt.xlabel("Year")
plt.ylabel("NO2 (ppm)")

Text(0, 0.5, 'NO2 (ppm)')

# O3
o3_df.plot(figsize=(40,25), color = plot_colors, linewidth=3.5, style = plot_line_styles)
plt.title("Ozone Trends from 2000 to 2016")
plt.xlabel("Year")
plt.ylabel("O3 (ppm)")

Text(0, 0.5, 'O3 (ppm)')

# CO
co_df.plot(figsize=(40,25), color = plot_colors, linewidth=3.5, style = plot_line_styles)
plt.title("Carbon Monoxide Trends from 2000 to 2016")
plt.xlabel("Year")
plt.ylabel("CO (ppm)")

Text(0, 0.5, 'CO (ppm)')

# Total Mean
total_df.plot(figsize=(40,25), color = plot_colors, linewidth=3.5, style = plot_line_styles)
plt.title("Total Gas Trends from 2000 to 2016")
plt.xlabel("Year")
plt.ylabel("Total Mean (ppm)")

Text(0, 0.5, 'Total Mean (ppm)')

Spread of Data¶

One would now like to know which of the states have most variance (spread/standard deviation) in the values from 2000 to 2016. In order to accomplish this, an ordered dictionary has been created with keys as state names and corresponding values as the standard deviation of their values.

def get_std_dict(df):
    """
    Returns an ordered dictionary with keys as states and values as standard deviation of gas emissions from
    2000 to 2016.
    
    Args:
        df (dataframe): corresponds to one the gas dataframes computed using get_df(gas_name) function.
    
    Returns:
        sorted_dict (dict): dictionary with keys as states and values as standard devations of gas emissions.
                            It is sorted in ascending order of values.
    """
    std_dict = {col:df[col].std() for col in df.columns}
    sorted_dict = sorted(std_dict.items(), key=lambda kv: kv[1])
    return sorted_dict

get_std_dict(total_df)

[('Pennsylvania', 0.062188213519883906),
 ('Oklahoma', 0.08541836974412767),
 ('New York', 0.09177862518857477),
 ('Texas', 0.095434416761091),
 ('California', 0.1116326201193114),
 ('North Carolina', 0.11648243132316333),
 ('Colorado', 0.11762492000856321),
 ('Kansas', 0.11764702515294129),
 ('Florida', 0.13268996670718286),
 ('Arizona', 0.1393261490093493),
 ('Illinois', 0.1433897740712859),
 ('Indiana', 0.19897774233960153),
 ('Louisiana', 0.22033401784133239),
 ('District Of Columbia', 0.3461618455293802)]